%% Step 4

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% AUTHOR: Produced for by JDG/IM based on initial code for cobb-douglas
% production in "Firm Model" directory.
%
% DATE: February 2019.
%
% Purpose: Computes the main algorithm to work out wages, prices, 
% endogenous labour for calibrations with shocks.
%
% Output:  converged expenditure hat, wage hat and prices, 
%         which are calculated via EquilibriumCES.m, which calls up the 
%         following functions:
%           1) PriceLoopCES.m
%           2) TradeShareCES.m
%           3) ExpenditureCES.m
%           4) LabourMarketConditionCES.m
%
% Last Updated: 07/09/2019 by JDG/IM
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

global alphaT rhoT lambdaT etaT shockT rho sigma lambda eta ...
    tol maxit Pi_hat_n D_hat_n Xi_hat_mnjf Xi_hat_mnj a_hat_f a_hat ...
    varphi sL_n sPi_n sD_n N J FI_J france oneminus_alpha_vector ...
    labour_share_PWT countrysecD firmsecD_sorted start start_sorted ...
    sum_by_country_dummy ISO alpha_WIOD_j alpha_WIOD_francej sector_algorithm 


%% Algorithm

tol      = 1E-5; %Tolerance allowed until price or wage converged
maxit    = 1E+10; %Maximum iterations to loop over

%% Initialise price hat, wage hat and expenditure hat guesses:
 
% P_hat: Initialise the Price hat equation 

P_hat_mnj0=ones(J*N,N);

% Initialisation of wage shock

w_hat0= ones(1,N);

% Initial guess for the change in expenditures 

X_hat_mnj0=ones(N*J,N);

%% Load in the updated hat data, taking into account the wedges

% Load in the data after completing the convergence the first time, 
% therefore calculating the new hat variables when the wedge is equal to 0. 
% To do counterfactual exercises we start from this point, therefore
% allowing our variables at time t to be equal to the hat variables
% calculated previously. 

%eval([strcat('load MAT/rho',rhoT,'eta',etaT,'lambda',lambdaT,'/Step_4_setup_alpha',num2str(alphaT),...
 %    '_rho',rhoT,'_approx.mat')]);  
% Function to solve the equilibrium 

[P_hat_n,P_hat_nj,P_hat_mnj,w_hat,X_hat_mnj,pi_mnjf1,pi_mnj1,...
    pi_c_mnj1,pi_c_nj1,pi_M1,pi_l1,pi_M_f1,pi_l_f1,PC_hat_n1,...
    final_goods1,intermediates1]= ...
    EquilibriumCES_fun_approx(P_hat_mnj0,w_hat0,pi_mnj,pi_M,pi_M_f,...
    pi_l,pi_l_f,pi_mnjf,pi_c_mnj,X_hat_mnj0,pi_c_nj,PC_n,zeta_mnj,X_mnj);